{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 15\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The coffee cooling problem\n",
    "\n",
    "I'll use a `State` object to store the initial temperature.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>T</th>\n",
       "      <td>90</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "T    90\n",
       "dtype: int64"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "init = State(T=90)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a `System` object to contain the system parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>T    90\n",
       "dtype: int64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>volume</th>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>r</th>\n",
       "      <td>0.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T_env</th>\n",
       "      <td>22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_0</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "init      T    90\n",
       "dtype: int64\n",
       "volume                     300\n",
       "r                         0.01\n",
       "T_env                       22\n",
       "t_0                          0\n",
       "t_end                       30\n",
       "dt                           1\n",
       "dtype: object"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coffee = System(init=init,\n",
    "                volume=300,\n",
    "                r=0.01,\n",
    "                T_env=22,\n",
    "                t_0=0,\n",
    "                t_end=30,\n",
    "                dt=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The update function implements Newton's law of cooling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_func(state, t, system):\n",
    "    \"\"\"Update the thermal transfer model.\n",
    "    \n",
    "    state: State (temp)\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: State (temp)\n",
    "    \"\"\"\n",
    "    r, T_env, dt = system.r, system.T_env, system.dt\n",
    "    \n",
    "    T = state.T\n",
    "    T += -r * (T - T_env) * dt\n",
    "    \n",
    "    return State(T=T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how it works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>T</th>\n",
       "      <td>89.32</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "T    89.32\n",
       "dtype: float64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "update_func(init, 0, coffee)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a version of `run_simulation` that uses `linrange` to make an array of time steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Runs a simulation of the system.\n",
    "    \n",
    "    Add a TimeFrame to the System: results\n",
    "    \n",
    "    system: System object\n",
    "    update_func: function that updates state\n",
    "    \"\"\"\n",
    "    init = system.init\n",
    "    t_0, t_end, dt = system.t_0, system.t_end, system.dt\n",
    "    \n",
    "    frame = TimeFrame(columns=init.index)\n",
    "    frame.row[t_0] = init\n",
    "    ts = linrange(t_0, t_end, dt)\n",
    "    \n",
    "    for t in ts:\n",
    "        frame.row[t+dt] = update_func(frame.row[t], t, system)\n",
    "    \n",
    "    return frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's how it works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>T</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>89.32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>88.6468</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>87.9803</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>87.3205</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>86.6673</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>86.0207</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>85.3804</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>84.7466</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>84.1192</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>83.498</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>82.883</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>82.2742</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>81.6714</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>81.0747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>80.484</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>79.8991</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>79.3201</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>78.7469</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>78.1795</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>77.6177</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>77.0615</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>76.5109</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>75.9658</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>75.4261</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>74.8919</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>74.3629</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>73.8393</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>73.3209</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>72.8077</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>72.2996</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          T\n",
       "0        90\n",
       "1     89.32\n",
       "2   88.6468\n",
       "3   87.9803\n",
       "4   87.3205\n",
       "5   86.6673\n",
       "6   86.0207\n",
       "7   85.3804\n",
       "8   84.7466\n",
       "9   84.1192\n",
       "10   83.498\n",
       "11   82.883\n",
       "12  82.2742\n",
       "13  81.6714\n",
       "14  81.0747\n",
       "15   80.484\n",
       "16  79.8991\n",
       "17  79.3201\n",
       "18  78.7469\n",
       "19  78.1795\n",
       "20  77.6177\n",
       "21  77.0615\n",
       "22  76.5109\n",
       "23  75.9658\n",
       "24  75.4261\n",
       "25  74.8919\n",
       "26  74.3629\n",
       "27  73.8393\n",
       "28  73.3209\n",
       "29  72.8077\n",
       "30  72.2996"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = run_simulation(coffee, update_func)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the results look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3hVVdbA4d9NJRA6hA4BhEWH0BGQEkWxix1UdBxBaVIVccb6KYgU6UUdwTYOiIplUDQQpApKb4saICH03knI98e5wZhJwgWS3Nxkvc+TJ7nnnrPPuhzNyt5nnb1dSUlJGGOMMTmNn7cDMMYYY9JiCcoYY0yOZAnKGGNMjmQJyhhjTI4U4O0AspuIBANNgHgg0cvhGGNMXucPlAFWqOr5lG/kuQSFk5wWejsIY4wxf9EaWJRyQ15MUPEAn332GaVLl/Z2LMYYk6ft27ePLl26gPt3c0p5MUElApQuXZry5ct7OxZjjDGO/7nlYkUSxhhjciRLUMYYY3IkS1DGGGNyJEtQxhhjcqRsLZIQkfbAu0B1YCswWFXnikgQMB54AOdG2ShVHZpOGy7gTaAbEAR8BAxS1YRs+AjGmDzk0qVLxMbGcvr0aW+H4rMCAwMJCwujUKFCV31stiUoEQkHvgMGA5Nxnkf6VkTaAI8BAlQFCgM/ikicqn6cRlPdgE5AQ+A88DUwBHgjqz+DMSZvOXToEC6XCxHBz88GnK5WUlISZ8+eJS4uDuCqk1R2/ot3BDar6jhVvaiqS4CZwJNAV+AtVT2qqjHACKB7Ou10Bd5T1VhVPQi8lsG+mW7vwVOcOnsxu05njPGiY8eOUapUKUtO18jlcpE/f37KlSvHgQMHrvr47PxX9wNS95MTcYb7ygAbU2zfDNRNp51aaexbVkSKZVKc6Tp68hzPvhNF96G/sGLjvqw+nTHGyxITEwkMDPR2GD4vJCSEixev/g/77ExQc4CGItJFRAJFpDnwCFDQ/f6ZFPueAfKn005oGvuSwf6ZplCBYBrVKMWJ0xd448Pf+PDb9VxMuJTVpzXGeJHL5fJ2CD7vWv8Nsy1BqeoO4B6gH7AP577RdCA5rYak2D0/cCqdpk6nsS8Z7J9p/P1c/PNvzeh6Ry38/Vx8s2A7L4z7lb2HsvzUxhiT52RbghKRgsBBVW2sqsVV9W6gLBCFk7Akxe41+OswXkob09g3XlWPZUHY/8PPz8UD7asxrFcrworlZ1vscfqOiib6jz3ZcXpjjMk0Q4cOJSIigltvvZUDBw5w3333ERERwdixY70dGpC9ZebFgd/cpeYrcHpTHXB6VCWAV0VkLc4Q3kBgTDrtfAIMFJEonN7Ua+5t2apGpWKM6d+W8TNXs3jNXkZ+vpJVWw7ybKd6hATnxSkOjTG+5tNPP+XDDz+kefPmfPvtt5w8eZIVK1YQEJAzfodl5xBfDPA0TjI5CgwAblfVfcArwHpgA07ymoVTig6AiJwSkS7ul5Nxqv+W4DxLtdF9fLYLDQnkxccb0+vB+gQF+jPv9z30Gx3Njrjj3gjHGJOHrVq1iocffvhyj2ju3LkkJCQwduxY2rRpQ7NmzXj22WeJjY0lMTGRiIgIEhIS6N69O+3bt+fll18mLi6OJk2asHPnTvbt20fPnj1p1qwZN998M9OmTbt8rsTERCZPnkxkZCTNmjXj+eef58iRI5n+mVxJSUmZ3mhO5n4ea2dUVFSmzma+e98Jhn/yO7v2nSTA34+/3VWbO1tVthusxviwTZs2UbNmzcuvX/9gGb9v2p8t525csxSv/r25R/seOXKEDh060K9fPx5++GFWrFhB9+7dueOOO1i3bh2TJk0iLCyMd999l2XLlvH1118TGBiIiPDNN99Qs2ZNvvrqK6ZPn87s2bNJTEzk/vvvp3HjxgwaNIj4+Hi6d+/Oc889x7333stHH33EzJkzmTx5MmFhYbzzzjvs2LGD6dOnpxtj6n/LZLGxsURGRgJUdndkLrPi/kxSsXQhRvZtQ8cW4SQkXmLqN+t466PlnDh9wduhGWNyufnz51OqVCm6dOlCQEAALVq04PPPP+eXX36hR48eVKhQgeDgYF544QX27t3L2rVrM2xv/fr17N69m8GDBxMcHEx4eDhPPfUUX3zxBQAzZsygV69eVKxYkXz58jFo0CBWrFhBTExMpn6unDHQmEsEB/rT44H61K9eknEzVvPbhn30GTmfAZ0bUfeGEt4OzxhznTzt0WS3w4cPU6ZMmb9sq1OnDmfPnqVs2bKXtwUFBREWFsa+fRk/xxkXF8fZs2dp3vzPz3vp0iWKFCkCwN69e3n55Zd55ZU/764EBAQQFxdHeHh4Jnwid5uZ1pK5rGW9slQrX4QRn/3BppgjvDx5MQ9FVufRDoK/v3VajTGZKywsjP37/zr0OG3aNJKSkoiLi6NBgwYAXLhwgf3791O8ePErtle8eHEWLfpzBfYjR45w7ty5y++/8sortG7d+vL7W7ZsydTkBDbEl2XCiuVnaI+WPHxLdQD+88sWBk9YxP4jZ65wpDHGXJ02bdqwf/9+Zs6cSWJiIkuXLmXs2LH06NGDSZMmsWfPHs6fP8/w4cMpWrQoDRs2zLC9evXqERoaysSJE7lw4QJHjhyhR48el8vP7733XiZMmEB8fDyJiYlMnTqVLl26XE5gmcV6UFnI39+Px26rSf0bSjLy8z/YvOsoz4+cT88HGtA6opy3wzPG5BJFixZl6tSpDB06lGHDhlGqVClGjhxJq1atuHjxIk888QTHjx+nYcOGfPTRRwQFBWXYXlBQEFOnTuXtt9+mdevWuFwubr75ZoYMGQJAt27dSEhIoEuXLhw7dozq1avz4YcfXtOM5RmxKr5scuL0BcbNWMWy9c7Y7y1NK9Lt3rrks2emjMmx0qs8M1fPqvhysEIFghjyZFOeu78eQQF+/Lx8N31HL2B7bLZMgGGMMT7HElQ2crlc3H5jZUb1bUPF0gWJO3iKgWMXMvvX7eS1nqwxxlyJJSgvqFSmEKP6tqHjjc4zUx/MXs/rHyzj2Mnz3g7NGGNyDEtQXhIc6E+P++sz5MmmhIYE8sfmA/QeOZ+Vm69+US9jjMmNLEF5WYu6ZRg3sB11q5bg2MnzvPr+Uvc6U4neDs0YAzb8ngmu9d/QElQOUKJICG8+eyNP3F4TP/c6UwPHLGTP/pPeDs2YPM3f3/+aVoI1f3X27NlrWpnYElQO4e/n4sHI6gzv1YrSxfOzY+9x+o5ewI9LY+wvOGO8pEiRIuzfv59Ll2zl7GuRlJTEmTNniIuLIyws7KqPt4dwchhxrzM15et1zPt9DxO+XMNKPUCvBxtQqEDGD9cZYzJXiRIliI2NRVW9HYrPCgwMpFSpUtf0EK8lqBwof75A+j3akIYSxsRZa1i6Lp4tu4/Sv3ND6t1Q0tvhGZNn+Pn5UbFiRW+HkWfZEF8O1qZhecYOaEfN8GIcPn6Of0xewvQfNpKQaMMNxpjczxJUDlfKPels5w6CC/hy3lZeGLeQvQdPeTs0Y4zJUpagfIC/vx+P3lqDoT1bEVY0hK17jvH8qGjm/rbLCiiMMbmWJSgfUqtyccYMaMdNEeU4dyGRcTNWM3T6Clu11xiTK2VrkYSINAfGAgIcBIap6gciknq8KgAIBsqp6t402tkNFAeSuw9xqipZF3nOERoSyKDHGtOkZikmzlrL0nXx6K6j9H+0IfWrWwGFMSb3yLYEJSJ+wGxggKp+KiJNgIUiskJVQ1PsFwDMB6LTSU4lgHJAIVU9nU3h5zhtG1WgRngxRn2+kk0xR/jHlCXc1/YGHu9Yg8AAf2+HZ4wx1y07h/iKAmGAS0RcOL2fBCD1+NSLQCDwajrtNAK25uXklKx08QIM7dGSLrfVwM/PxdfR22wGCmNMrpFtCUpVDwPjgenARWAFMERVNyXvIyJlgSHAs6qaXi11Q8BPRJaLyEER+UlE8uyKYv7+fjxyi/BOyhkoRkXz3yU7rYDCGOPTsi1BuYf4zgGdgRCgLfCqiHRIsVs/4EdVXZ1BU4nAcqATUAlYBcwRkfxZEbevqOGegSKySQUuJFxi0qy1vPmv32wJD2OMz8rOIb5OQEtV/UJVL6rqAuBDoDuAiPgDXYHJGTWiqsNV9TFVjVXVM8BLQDGcob88LX++QPo+0pAXHm9MgZBAVmzcT+8R8/l9035vh2aMMVctOxNUBZzKvJQScIb7AG50f4/KqBER6SsirVJs8scp9jiXGUHmBq0blGPcAPcSHqfO8/oHy5g0aw3nLiR4OzRjjPFYdiaouUBtEekmIi4RaQQ8A/zb/X5zYFkG956ShQNjRKSMiIQAo4CtwMositsnlSzqLOHx1J21CPB38d8lMfQdtYBte455OzRjjPFIdhZJbMAZ5usOHAM+Bwar6mz3LuHA/5SVA4jIKRHp4n45GFiGc+/pAFAFuEtVbYW/VPz9XHRqV42Rz7ehQqmCxB08xcCxvzIzaguJl6yAwhiTs7nyWqWXiIQDO6Oioihfvry3w8k25y8mMu37DXy/aCcAtSoXo3/nRpQqlqdrS4wxXhYbG0tkZCRAZVWNSfmeTXWURwQH+tP9vnq89kxzihYMZuPOI/QeMZ95v++xcnRjTI5kCSqPaVSjFOMGtqN5ndKcPZ/A6H+vZPgnv3PyjM3nZ4zJWSxB5UGFQ4MZ8mRT+jzUgHxB/ixas5feI+azZstBb4dmjDGXWYLKo1wuF7c0q8TYAe2QSkWdBRGnLOGD2eu5cNHqTYwx3mcJKo8rU6IA7/RsRedbnfn8Zv+6nX7vLWBH3HFvh2aMyeMsQRlnQcQOwru9W1O2RAF27zvJgDEL+HLeVitHN8Z4jSUoc1n1ikUZ078tt98YTkJiEtN/2MjLkxaz/8gZb4dmjMmDLEGZv8gXHMBz99fn1b875egbdhym94j5/LJ8t5WjG2OylSUok6bGNZ1y9BZ1y3D2fAJj/rOKodNXcPyUzY5ujMkelqBMugqHBvNS1yb0fSSCkOAAlq6Lp5fNjm6MySaWoEyGXC4XkU0qMm5gO2pXKc6xk87s6BO/XMO58zY7ujEm61iCMh4pVSw/bz3X8vLs6HOWxtBnVDSbY454OzRjTC5lCcp4LHl29FF92xBephDxh07z4viFfPzfjVxMuNIqKcYYc3UsQZmrVrlsYUb1vYn7291AEjAzaisDx/zKrvgT3g7NGJOLWIIy1yQwwJ8n76zN0B6tKFUsPzv2Hqfv6AV8NX+bPdxrjMkUlqDMdaldpThjB7Tl1uaVSEi8xEffb+DlSYvZd/i0t0Mzxvg4S1DmuuXPF0ivBxvwytPNKOJ+uLfPyPn8/Nsue7jXGHPNAjzZSUQ6Ah2BxkAYkAjsA1YA36tqdFYFaHxHk1qlGT+wHRNnrWHJ2njGzljNbxv20fPB+hQtmM/b4RljfEyGCUpEHgVeA4oBvwBzgMOAP1ACqA/MEJGDwFuq+nmWRmtyvMKhwQx+ognRK2OZ8tVaftuwj00xR+jxQH1a1ivr7fCMMT4k3QQlIj8Bp4FngQWqmmYdsYj44fSunhORp1T1lgzabA6MBQQ4CAxT1Q9EJBg4CaRc1nWJqnZIow0X8CbQDQgCPgIGqao9NZpDuFwu2jWqQJ0qJRj7n1Ws3nqQYdNX0CaiPM92qkto/iBvh2iM8QEZ9aBeUdXfrtSAO3H9APzgTkBpciey2cAAVf1URJoAC0VkBRAIHFHV0h7E3A3oBDQEzgNfA0OANzw41mSjkkVDeL1bC+Ys2cm/vt/IglWxrNt+iD4PN6BRjVLeDs8Yk8OlWyThSXJK45hlGbxdFOf+lcvdC0oCEnB6TY2A1R6epivwnqrGqupBnCHI7lcbq8kefn4u7mhVhXED2lKjUlGOnDjHa+8vY8KXazhrUyUZYzKQYRWfiFQVkU9FpFKq7dNE5IvU2zOiqoeB8cB04CJOgcUQVd2E0xsKE5G1IrJfRGaKSLl0mqoFbEzxejNQVkSKeRqLyX5lS4YyrFdrut5RiwB/P35cGkPvEfPZsOOwt0MzxuRQ6SYoEakGLMUphCiQ6u35QE3gNxGp7MmJ3EN854DOQAjQFnhVRDrg3OtaDETi3J86izN0l5ZQIOUKesk/5/ckDuM9/n4uHmhfjdH92lC5bCH2HznDSxMX8eG367lwMdHb4RljcpiMelBv4iSoCFVN2WNBVacDTYB1eH7vpxPQUlW/UNWLqroA+BDorqr9VbW3qh5U1WNAf6CJiFRIo53TOAkuWXJiOuVhHMbLwssUYuTzbXjo5uq4gG8WbKfv6Gi27jnq7dCMMTlIRgnqJuCN9KrjVPUCTnJq5+G5KgDBqbYlABdF5A0RqZlie3KZ17k02tmI08tKVgOIdyc24yMCA/x4vGNNhvduTbmSoezZf4qBYxfy2Y+bbeJZYwyQcYIqCFzpT9pYoLCH55oL1BaRbiLiEpFGwDPAv4F6wEgRKSIiRYAxwA/uIojUPgEGikglESmBUyTxiYcxmBxGKhXjvf5tuLt1FS5dSuKLn5WBY34lxiaeNSbPyyhBbQLSLRt3awHEeHIiVd2AM8zXHTgGfA4MVtXZwNM4yXCbu70LwOPJx4rIKRHp4n45GZgJLAG24vSoXvEkBpMz5QsK4Jl76/J2j5aXJ57tNzqamVFbSEy03pQxeZUrvbnSRKQbzvNFbVU1Jo33qwDRwDhVfTcLY8xUIhIO7IyKiqJ8+fLeDsekcvZ8Ah99t4E5S2MAqF6xCH0faUiFUgW9GpcxJmvExsYSGRkJUDl1rkn3QV1VnSoibYGNIvIlTln4cZznmZrg9IbmAKOyJmyTF4UEB9Djgfq0qFuGsf9ZxZbdx+g7KprHb6/JXa2r4u/n8naIxphskuFzUKraGfg7UBZ4FfgXTq+qEPCYqt6vqlYfbDJdhIQxflB7bm5SkQsJl/jw2w0MmbiI+EO2jIcxecUVZzN3TwBrk8CabFcgJJDnH4mgRb0yjJ+xmo07j9B75HyeurM2HVuE42e9KWNytYwe1H1ZRELSez+N/QuIiBUrmEzXtFZpxg9qT5uI8py/kMjkr9byytQlHDhy5soHG2N8VkZDfOeB9SLytog0Tm8nEYkQkZE4VX9pPbdkzHUrVCCIgY81YnDXJhQODWLN1kP0GjGPOUtjbFFEY3KpjIokRojITOBFIFpEzuKUdB8CXEBJoA7OQ7UfAzelVe1nTGZqWa8sdaoUZ9JXa1m8Zi8Tv1zDkrV76f1QA8KK2mxXxuQm6ZaZpyQiBYH2OLOOlwIu4ayo+wcQpapnszLIzGRl5rnHwtVxTJq1lpNnLhASHMDTd9emQ7NKuFx2b8oYX3FNZeYpqepJnLWcZmd6dMZco9YNylG3agkmfeUsMT9+5hoWr9lL74ciKFnU49unxpgcKsMyc2NyuiIFnSXmBz3WiIL5A1m15SC9Rsxj7m+77N6UMT7OEpTxeS6Xi5siyjPhhfa0qFuGM+cSGDdjNa99sIxDx3xm9NkYk4olKJNrFC2Yj5e6NmFAl0aEhgSycvMBer47j5+tN2WMT7IEZXIVl8tF24ZOb6pZ7dKcOZfA2Bmree39ZRw8ar0pY3yJR0USAO5nofoD1YF7gYeBHaqa3sq3xnhNsUL5ePmppkSvjGXq1+tYqU5vyir9jPEdHvWgRORWYAHOUuy1cZ59yg/8R0SezLLojLkOLpeLdo0qMPGF9jSvU5qz5xMYP3MNr0xdarNQGOMDPB3i+z9ggKo+jbMKLqr6JtAH50FeY3KsooXyMeTJpu5KvyBWuyv95izZyaVLdm/KmJzK0wRVC2dF3NTmAuGZFo0xWeTPSr923FivDGfPJzJx1lr+OWUJ+w7bDOnG5ESeJqg4ICKN7e2B3ZkXjjFZy6n0a8qLTzSmUIEg1m47RO8R8/lh0Q7rTRmTw3haJDEUmOpeRdcfuM09ZVBPoF8WxWZMlmlV35mFYvJXa1m0Zi+Tv17HorV76fNQBGVKFPB2eMYYPOxBqepHwBPAPTgzlg8D2gJdVXVqlkVnTBYqHBrMi080YXDXJhQJDWb99sP0GjGf2b9uJ9F6U8Z4nUc9KBF5AfhMVVtdz8lEpDkwFhDgIDBMVT8QkTBgDBCJM1P6HOB5VT2aTju7geJA8m+ROFWV64nN5F3JM6S//816FqyK5YPZ61m0Oo4+D0dQoVRBb4dnTJ7l6T2ol4F813MiEfHDmWx2rKoWBh4FxotIfeADnOrAykA1oCgwIZ12SgDlgDBVDXV/WXIy16VwaDADH2vEP55qSrFCwWzedZQ+I6OZGbWFxMRL3g7PmDzJ0wT1A9BbRIpex7mKAmGAS0RcOL2fBPfXJeB1VT2tqseA94H0emuNgK2qaqVXJtM1q1OGCS9EckvTiiQkXuLj/25iwNhf2bn3uLdDMybP8bRI4gbgEZwkdQLngd3LVLXslRpQ1cMiMh6YDnyEU2zxvKpuwJmZIqV7gVXpNNUQ8BOR5Tg9rpVAX1Xd5OFnMSZDoSGB9Hk4glYNyjFh5mq2xx6n3+gFPBhZnYdurk5ggM0QZkx28DRBTXF/XTP3EN85oDMwC7gR+EpENqvq3BT7DcRJUDem01QisBwYDBwBXgHmiEgtVbXpAUymaShhjBvYjo//u4kfFu/ki5+Vpev28vwjEVSrcD2DCcYYT3i0om5mEJEHgP6qemOKbcOBqqp6v4gEAuOAu4A7VHW1h+26gOPuYxZ6sH84tqKuuUrrtx9i7IzVxB86jZ8L7mt7A4/eWoPgQH9vh2aMT7vuFXVF5POM3lfVzh40UwEITrUtAbjoXlL+O6Ag0FRV4zKIpS/wu6oucm/yx/kc5zyIwZhrUqdqCcYOaMtnP27m21+3M2v+Npatj6f3QxHUrlLc2+EZkyt5OpiemOrLBVTFGYrzdCaJuUBtEekmIi4RaQQ8A/wb+MIdS+uMkpNbODBGRMqISAgwCtiKcy/KmCyTLyiAp++uw/DeralQqiBxB08zeMIiJs1aw5lzF70dnjG5jkc9KFV9PK3tIjIEJ1F50sYGEekEvAm8C+zDuY+0E7gdOA8cELlcMX5MVcu7z3MK6K6qn7mPGYlTRFEAZ5b1u1Q10ZM4jLleUqkYY/q3YcYvW5kZtYX/Lolh+cb99HygPo1rlvJ2eMbkGtd1D8p9P2etqhbKtIiymN2DMplp597jjJ2xmm17jgHQtlF5/n53HQqHph7NNsakJaN7UNdbL/skcPI62zDGZ1UuW5gRvVvzt7tqExTgR/QfsfR8dx4LV8fZMvPGXCdPiyTi+XNaoWSh7q8BmR2UMb7E39+P+9reQLM6pRk3YzXrtx9m+Ce/s2BlaZ67vx7FC4d4O0RjfJKnz0H9I9XrJOACsEJVNXNDMsY3lS0RylvPtmTub7v46PsN/LZhH+u3H+Kpu+rQoVlFW2bemKvkaYIqBYxW1b/MICEihURkuKq+kPmhGeN7/Pxc3NYinMY1SzFx1hpWbNzP+Jmr+XVVLD0frE/ZEqHeDtEYn5FughKRMjjPJYFTeTdfRA6n2i0C6A1YgjImhRJFQvjn35rx66o4pn6zzlkY8d35dL61Bve2qYq/v02XZMyVZNSDuhGYyZ/3nhaneC8J51kogGmZH5Yxvs/lctGmYXkaVC/JB9+uJ/qPWKb9sJFfV8fR56EGVC1fxNshGpOjpftnnKrOwpkkVnCSUQucpTCqAdXd75VW1b9lQ5zG+KzCocEM6NyI155pTljREHbEHaf/mF+Z9v0Gzl+0x/eMSU+G96BUdQeAiASm9yCs+z17jN6YK2hUoxTjB7Xn0zmb+G7RDmbN38aSdfH0frABdW8o4e3wjMlxPC2SKCYig4FaOHPfgdOrCgbqAMWyIDZjcp2Q4ACeubcurSPKMW7GanbvO8mQSYvp0KwST91Vm9CQQG+HaEyO4emd2qnAQzjTErXDmfvuEs6igm9lTWjG5F41KhXjvX5t6XJbDQL8/Zj72y56Do9i6bq93g7NmBzD0wTVHuiqqj2ADcC/VPVWYDjpr9tkjMlAYIAfj9wijOnfhprhxThy4jxvT1vB29OWc/j42Ss3YEwu52mCCsbpNQFswll2HeBfQMvMDsqYvKRi6UIM69mK7vfVJSTYn6Xr4ukxfB5zluzk0iWbLsnkXZ4mqC38mYg2Ac3cP+d3fxljroOfn4s7W1VhwqBImtYqzZlzCUyctZaXJi5iz36b7tLkTZ4WSYwApomIP/AfYI2IJAHNgSuuYmuM8UzJoiH8429NWbx2L1O+XsfGnUfoMzKahyKr8UBkNQIDbAVfk3d41INS1Y+Bm4FN7rn37sGp3FsKPJ114RmT97hcLlrVL8ekF9pza/NKJCRe4vO5yvOjotm4M/VkLsbkXp7OZj4beFFVNwOo6k/AT1kZmDF5XWj+IHo92IA2DcszYeZq9uw/xYvjF9HxxnC63l6LAlaSbnI5T+9BtcKZvdwYk83qVi3B2AHtePjm6vj7uZizJIYew+dZSbrJ9Ty9BzUKmC4io4AdwF9qYFV1S2YHZoz5U1CgP491rEnrBuUYN3M1uusob09bQfM6pel+Xz1KFLE1p0zu42mCetP9PbmSL7n21eX+2e7cGpMNKpUpxDu9WjNnyU4+/u8mlq3fx5qth3i8Y01ub1kZfz9bc8rkHp4mqGpZGoUxxmP+7pL05nXKMPWbdSxdF8/Ub9YRvXIPvR5sQOWyhb0dojGZwqMEparbAUSkOs7s5lFACVXdfTUnE5HmwFh3GweBYar6gYgEAeOBB4BEYJSqDk2nDRdOj64bEAR8BAxS1YSricUYX1eiSAhDnmzK0nXxTPl6LVt2H6Pv6AXc16Yqj3QQ8gV5+venMTmTR0USIlLQXcm3GfgaKA2MEZFVIlLWwzb8gNnAWFUtDDwKjBeR+sDrOEmrKtAE6CoiT6TTVDegE9AQp2fXBBjiSQzG5EYt6pZh4gvtubNlZZKSkpg1fxu93p3Pys0HvB2aMdfF0yq+EUBhIJw/CyT6A6eA9zxsoygQBrjcvaAkIAGnOrAr8JaqHlXVGPf5uqfTTlfgPVWNVdWDwGsZ7GtMnpA/XyDdO9Xj3d6tCS9TiP1HzvDq+0sZ8ekfHDt53tvhGXNNPE1Qdwa/XjsAAB5vSURBVAIDUg7pqepOoA8Q6UkDqnoYZxhvOnARWIHT84kHygAbU+y+GaibTlO10ti3rIjYkh8mz5NKxRjdrw1P3lGLoEB/FqyK5bl3opj72y6b18/4HE8TVEHgTDrvefqwrx9wDugMhABtgVeBu927pGz/DOnP8Reaxr5ksL8xeUqAvx/3t6/GhEHtiKheklNnLzJuxmqGTFrM7n0nvB2eMR7zNEH9BLzinosPIElEigPvAj972EYnoKWqfqGqF1V1AfAhzpAdOEkrWX6c4cO0nE5jXzLY35g8qXTxArzerQUDujSiSGgwG3Yc5vlR0XwyZ5MtNW98gqcJqg9QGTiEkxB+BHYDJYC+HrZRAWfZjpQScKr59uEUSSSrwV+H8VLamMa+8ap6zMM4jMkzXC4XbRuWZ9KLyfP6JTHjly30fnc+q9SKKEzO5mmZeTzQXERuwbkHFICz7MaPqnrJw3PNBYaKSDfgfZwqvGeAv+Mku1dFZC3OEN5AYEw67XwCDBSRKJze1GvubcaYdCTP69e+cQUmfLmG3ftO8srUpbSJKM/T99SmaMF83g7RmP/haQ8q2Q4gBlBg/VUkJ1R1A84wX3fgGPA5MFhVZwOvAOtxVutdAcwCJicfKyKnRKSL++VkYCawBGcRxY3u440xV1CrcnHe69eWJ26vSVCAn7uIYh4/LYuxIgqT47iSkq78H6WIFAFm4FTsJeJMceSH80zUE6qaXgFFjiMi4cDOqKgoypcv7+1wjPGafYdPM2nWWla6h/pqhhej54P1qVS6kJcjM3lJbGwskZGRAJXdjxld5mkPagrO/aamQD6ce0ktcB6s9fQ5KGNMDlK6eAFee6Y5LzzWmCIFg9kUc4TnR0Yz/YeNnLtgE7MY7/N0LpSOQFtVXZli23IR6Y5TMNEt0yMzxmQ5l8tF64hyRNQI4+MfNvLjshi+nLeVX1fH8VynejSuWcrbIZo8zNMe1CGcmSBSC8ApVDDG+LDQkEB6PFCf4b1bU7lsIQ4cOcPrHyxj6PTlHDp29soNGJMFPO1BvQRMFZGXgEU45eENgdHAOPcksoCtDWWML6tRqRij+7bhu0U7+OzHzSxZG88qPUCX22pyZ8vK+PtfbV2VMdfO0wT1b/f3L/jrWlAAw4Ch2NpQxuQK/v5+3NvmBlrWK8f7s53lPD6YvZ55K/bQ88H6VK+Y1mCKMZnP1oMyxqSpZFFnOY/lG/cx5au17Nh7nIFjf+W2FuE8cXstQkMCvR2iyeWuaj0oY0ze07RWaepVLcEXPyvfLNjOnCUxLF0Xz9N316FNRDlcLlvF12QNTyd6rQeMAmrzv9MVoao2k7gxuVi+4ACevLM27RpVYOKsNWzceYSRn/3BL8t38WynepQPK+jtEE0u5OkQ33Scar1/4MxIbozJgyqVKcTQHq2IWrGbj77fwJqth+g9Ipr7293AgzdXJzjQbkGbzONpgqoONFHV9CZwNcbkEX5+Lm5pVommtUsz/YeN/Lx8N//5ZQvRK2Ppfl9dmtQq7e0QTS7hac3oYpzhPWOMAaBwaDB9Ho5geK8/V/F948PfeHvacg4c9ZnZz0wO5mkP6hlgsYh0xJkw9i+TxKrq25kdmDHGN9SsXIz3+rXhu0U7+fynTSxdF89KPUDnDsLdN1UlwJ6dMtfI0wT1T5xl2VsCEaneSwIsQRmThznPTlWlVf2yfPDtehav2ctH328k6vc9PNepHnWqlvB2iMYHeZqgHgE6uZfGMMaYNJUoEsLgJ5rwx+b9TPlqHbv3neSliYtp37gCT91ZmyIF/6cI2Jh0edr3PgrYFEbGGI80qlGK8YPa8WgHIcDfj3m/7+HZd6L4YfFOEm3dKeMhT3tQLwJjRWQAsB24mPJNVb2Q2YEZY3xbUKA/nW+tQduG5Zny9TpW6gEmf7WWX5bv4rn7bcokc2We9qBGAm2AVcAJ4GyqL2OMSVPZkqG89kxzBndtQonC+dgW60yZNH7mak6ctr9tTfo87UE9lqVRGGNyNZfLRct6ZWkoYfzHPWXST8t2sWRtPF3vqMUtTSvi52dTJpm/8nQuvqjkn0WkIHBKVW0g2RhzVULcUyZFNqnI5K/WsnbbIcbPXM3Py3fR4/76VClX2NshmhzE0x4UItIfGASUBKqLyD+B48BAVb3i+tAi0gVn6fiUQoAo4MY04goGyqnq3jTa2g0U58+lP+JUVTz9LMYY76pQqiD/9+yN/Loqjn99tx7ddZR+o6O5vWVlHrutJgVspnSD55PFDgD6AC8Ak92bZwOTgAvu7RlS1c+Az1K0GQHMBQap6poU2wOA+UB0OsmpBFAOKKSqtpqvMT7K5XLRpmF5mtQqxWc/beb7RTv5ftFOFq3Zy1N31qJdowo2U3oe52mRxDPAs6r6Ce5ZJFT1G6Ar0PlqTyoigTjJ6rWUycntRSAQeDWdwxsBWy05GZM75M8XyDP31OW9fm2oVbkYx06eZ/S/VzF4wiJ27j3u7fCMF3maoCoCm9PYHgNcy1IbPXGq/yam3CgiZYEhOMnwUloH4iw17yciy0XkoIj8JCI1ryEGY0wOUrlsYYb1bEW/RyMoEhrMxp1H6Dsqmilfr+XU2YtXbsDkOp4mqFXAQyleJ9/7edb9nsdEJAjnXtZraRRa9AN+VNXVGTSRCCwHOgGV3OefIyL5ryYOY0zO43K5aN+4IpMGR3JX6yoAfL9oJ88NiyJqxW4u2UO+eYqnRRIDgR9F5Cac4oXXRaQGUAvoeJXnvA1nmPCHlBtFxB9nyLBLRger6vBUx70E9MAZ+lt4lbEYY3Kg0JBAut1bl1uaOtV+G3ce4b0vVvHTMmeBRKv2yxs86kGp6lKcNaFW4SSWokA0UFNVrzYp3APMSGMIL7mSL4oMiEhfEWmVYpM/TqK1hRSNyWX+HPZrSJGCwWyKOUK/0dFM+cqG/fKCdHtQIvIKMEJVzwCo6n6cFXWvV3Oc2dHT2r4sg3tPycKBx0XkTuAY8A6wFViZCbEZY3IYZ9ivAs1ql+bznzbz/eKdfL94JwvXxPHkHbVo39ge8s2tMupBvQqEZsE5w4H/KR/PYDsicsr9HBXAYGAZTm/uAFAFuEtVEzM9UmNMjlEgJJBn7v2z2u/4qQuM+c9qXhi/kG17jnk7PJMFXElJad90FJFLQGlVPZC9IWUtEQkHdkZFRVG+fHlvh2OMuQZJSUlEr4zlo+82cPTkeVwuuLV5OI93rEmhAkHeDs9chdjYWCIjIwEqq2pMyveuVCRRXkTyXekEqrr72sMzxpir43K5aNfIGfb791zlu4U7+HFpDIvXxPF4x5p0aB6Ovw37+bwrJagVV3jfhVNy7p854RhjjOfy5wvk6bvr0KFZJaZ+vY7VWw8ycdZafly2i2fvq0fNytfymKbJKa6UoNoBh7MjEGOMuVYVShXkje4tWLIung+/Xc+OuOO8MH4h7RtX4Mk7alG00BUHgkwOlFGCSgI257Z7UMaY3Cl5SY9GNcL4ct5Wvpq/jXm/72Hpung631qDO1tVJsDf07kJTE6Q0dWyAVxjjM/JFxTAY7fVZMKg9jStVZqz5xP48Nv19Bk5n9Vb7O9tX5JRgpqOrZZrjPFRZUoU4J9PN+PVvzenTIkC7Nl/in9OWcrb05az/8gZb4dnPJDuEJ+qPpWdgRhjTFZoXLMU9auVYPavO/jPz8rSdfH8sWk/ndpV4/72N5AvyONl8Uw2swFZY0yuFxjgzwPtqzF5cCRtG5bnQsIlvvhZ6TF8HovX7CW950GNd1mCMsbkGcULhzCgSyOG9WxFlbKFOXj0LMM+XsE/Ji9hV/wJb4dnUrEEZYzJc2pXKc6ofm3ocX89CuYPZO22Q/QZFc3Ub9Zx6swFb4dn3CxBGWPyJH8/Fx1vrMyUl27mjpaVISmJ7xbuoPuwKH5cGkOirT3ldZagjDF5WsH8QTzbqR7v9W9L7SrFOXH6AhO+XEP/9xawYYfNU+BNlqCMMQZn7amhPVrywmONKVEkhB1xxxk8YRHvfvI7B4/aEzfeYPWVxhjj5nK5aB1Rjia1SzFr3ja+mr+VX1fHsWzDPh6MrMZ9bW8gONCmHs0u1oMyxphU8gUF0OW2Gkx6MZKW9cty4WIin/24mR7vRFlZejayBGWMMekIK5afwU804e3nWhJephAHUpSlx1hZepazBGWMMVdQ94YSvNevDc+lKEt/fuR8Js1aw4nTVpaeVSxBGWOMB/z9/bjdXZZ+Z8vK4HLx3yUxdB/6C98u3E5C4iVvh5jrWIIyxpirUDB/EN071WNs/7Y0qFaSU2cv8v43zmzpKzfbbOmZKduq+ESkCzAl1eYQIAq4CzgJpOwrL1HVDmm04wLeBLoBQcBHwCBVTciKuI0xJi2VyhTije4tWL5hHx9+t4E9+0/x6vtLaVyzFH+/pw7lSoZ6O0Sfl20JSlU/Az5Lfi0iEcBcYBBQFziiqqU9aKob0AloCJwHvgaGAG9kdszGGJMRl8tFszplaFgjjO8W7uCLn7fw+6b9rNID3NW6Cg/fIoSGBHo7TJ/llSE+EQnESVavqeoaoBGw2sPDuwLvqWqsqh4EXgO6Z0mgxhjjgcAAfzq1q8aUlyLp0KwSl5KS+GbBdroP/cWmTboO3roH1RNnMcSJ7tcNgTARWSsi+0VkpoiUS+fYWsDGFK83A2VFpFjWhWuMMVdWtGA+ej/UgNF92/xl2qR+o6NZt+2Qt8PzOdmeoEQkCGdY7zVVTf6z4jSwGIgEBCd5fZ1OE6FAyuUwk3/On/nRGmPM1atavghDe7TkxScaE1Y0hJ17TzBk0mLenrac+EOnvR2ez/DGVEe3AZeAH5I3qGr/lDuISH/goIhUUNU9qY4/jVNckSw5MZ3KgliNMeaauFwuWtUvR5Napfkmehsz521l6bp4Vmzcz92tq/DQzdUpYPenMuSNIb57gBmqevmhARF5Q0RqptgnyP39XBrHb8TpZSWrAcSr6rFMj9QYY65TcKA/D98iTBkcSfvGFUhIvMRX0dvoPszuT12JNxJUc5zhvJTqASNFpIiIFAHGAD+4iyBS+wQYKCKVRKQETpHEJ1kZsDHGXK/ihUPo92hDRvW9iZrhxTh+yrk/1XdUNGu2pPWrzngjQYUDe1Ntexo4CmwDYnCeh3o8+U0ROeV+jgpgMjATWAJsxelRvZKlERtjTCapVqEo7/RqxQuPO/enYuJP8I8pS/i/f/3G3oN2pyIlV16blVdEwoGdUVFRlC9f3tvhGGPysPMXE/n21+3MjNrC2fOJBPi7uLNV3np+KjY2lsjISIDKqhqT8j2b6sgYY7wkONCfByOrM3nwzdzStCKJl5znp7q9/Qs/LNpBYh6f388SlDHGeFmxQvno83AEo9zPT508c4HJX6+j98j5/L5pf55df8oSlDHG5BA3uJ+fGvJkE8oUL8Ce/ad4/YNlvDJ1aZ5cf8oSlDHG5CAul4sWdcsy4YV2PH13bQrkC2D1loM8P3I+42eu5ujJtJ6+yZ0sQRljTA4UGODPvW1u+Mv6Uz8t20X3oVHMjNrChYuJ3g4xy1mCMsaYHKxwaDDdO9Vj/MB2NK5ZirPnE/j4v5t47p0oFq6Ky9X3pyxBGWOMD6hQqiCv/r05b3RrQXiZQhw4epbhn/7OC+MWsjnmiLfDyxKWoIwxxodESBjv9W9LrwcbUKRgMJt3HWXQuIW88/EK9h3OXRPRemOyWGOMMdfB38/Frc0r0bpBWWbN38Y30dtYtGYvy9bv4y73RLS54UFf60EZY4yPyp8vkMc71mTy4Jtp16g8CYmX+Dp6G93e/oXvFu4gwccf9LUEZYwxPq5k0RD6d27E6L5tqFPVedB36jfr6PXuPJatj/fZQgpLUMYYk0vcUKEIbz/Xkpefakq5kgWIO3iatz5azpBJi9m2x/dWJLIEZYwxuYjL5aJ5nTKMH9SebvfWpWD+INZvP0y/9xYw8vM/OHD0zJUbySEsQRljTC4U4O/HXa2rMHXIzXRqewMB/n5E/xHLs8OimPb9Bk6fvejtEK/IEpQxxuRioSGBPHVXbSa92J6bGpTjYsIlZs3fRrehTiHFxYScW0hhCcoYY/KA0sULMOjxxox8/iZqVynOidNOIUXPd+exZO3eHFlIYQnKGGPykOoVizK0x5+FFPGHTjN0+gpeHL+Izbty1owU9qCuMcbkMcmFFI1rluKnZbv499zNbIo5wqCxC2lZvyxdb69FmRIFvB2mJShjjMmrAvz9uKNlZdo1Ks+X87Yye8F2Fq/Zy2/r47m9ZWUevlkoVCDIe/Fl14lEpAswJdXmECAKeAwYA0QCLmAO8LyqHk2nrd1AcSB50DROVSUr4jbGmNwuf75Anri9Fh1bVObTHzcx/489fPvrDqKW7+aByOrc1boKwYH+2R6Xy1s3xkQkApgL3Ay8CRwHngUCgU+Ak6raOY3jSgD7gUKqetUzI4pIOLAzKiqK8uXLX/sHMMaYXGpH3HGmfb+BVVsOAlCicD4e61iTto0q4O/nytRzxcbGEhkZCVBZVWNSvueVIgkRCQQ+A14D1gGXgNdV9bSqHgPeB1qlc3gjYOu1JCdjjDFXVqVcYd7ofiOvd2tB5bKFOHT8HO99sYp+o6NZqQeyLQ5v3YPqCZwFJqpqEnBvqvfvBValc2xDwE9ElgOVgZVAX1XdlFXBGmNMXtRQwqhfrSQLVu7hkzmb2bn3BK9OXUqD6iV56s7aVClXOEvPn+09KBEJAgYBr7mTU+r3B+IkqBfTaSIRWA50AirhJLI5IpI/ayI2xpi8y9/PRfvGFZk8OJIn76hFgXwBrN5ykL6joxmVxVMneaMHdRvOkN4PKTe6h/3GAXcB7VV1c1oHq+rwVMe9BPTAGfpbmBUBG2NMXhcc6M/97atxS7NKzPhlCz8s3sH8P2JZvDae0X1vomLpQpl+Tm/cg7oHmKGql+fXEJGCwM9AE6Cpqq5O72AR6SsiKe9P+eMk2nNZFK8xxhi3QgWC+Ps9dZj0YiQ3RZSjSGgQfplcOJHMGz2o5sA/U237AidZtlbVK/UXw4HHReRO4BjwDrAV516UMcaYbFC6eAEGPdY4S8/hjR5UOLA3+YWI1ANuB5oCB0TklPsrNsU+p9zPUQEMBpbh3Hs6AFQB7lLVxGyK3xhjTDbI9h6UqhZI9XotzsO5GR0TmuLnczhVgD2zJEBjjDE5gk0Wa4wxJkeyBGWMMSZHsgRljDEmR7IEZYwxJkeyBGWMMSZHyovrQfkD7Nu3z9txGGNMnpfid/H/rOeRFxNUGYAuXbpcaT9jjDHZpwywPeWGvJigVgCtgXiciWeNMcZ4jz9OclqR+g2vLVhojDHGZMSKJIwxxuRIlqCMMcbkSJagjDHG5EiWoIwxxuRIlqCMMcbkSJagjDHG5EiWoIwxxuRIlqCMMcbkSHlxJolrJiL1gclAPWAH8DdV/Z+nn32NiPwNmAKcT7G5p6pO91JI10VEmgLfq2qY+3UQMB54AGf2kFGqOtSLIV6TND5XMHASuJBityWq2sEb8V0tEbkFGAZUAw4A76rqFF+/Xhl8Lp++XgAicifwNlAZ57MNz8prZgnKQ+4LMBt4D7gJuB+YKyKVVPWEV4O7fg2Bkao62NuBXA8RcQFPAyNSvfU6IEBVoDDwo4jEqerH2RziNcngc9UFjqhq6eyP6vqISAVgFtAV5/+rRsBPIhIDtMVHr9cVPtdhfPR6AYhIGeBL4D5VnSMiDYHFIrICeJAsuGY2xOe5tkCgqr6nqhdV9QtgA/Cwd8PKFI2A1d4OIhO8DjwH/F+q7V2Bt1T1qKrG4Pyi757NsV2P9D6XL1+3cOBzVf1aVS+5RyKigZb49vUKJ/3P5cvXC1WNB0q6k5MfUBxIwOkVZsk1sx6U52oBm1Jt24zzV6zPEhF/nCHLx0VkFHAG+AB4R1V9baLGyar6ioi0Td4gIkVwJqLcmGI/X7tu//O53BoCYSKyFigF/Ar0VdW47A7waqnqQmBh8msRKYYzifMn+PD1usLnug0fvV7JVPWkiOQHjuPkj3eAg2TRNbMelOdCcX55p3QGyO+FWDJTSeB3YDrOuPIDOH+tP+fNoK6Fqu5NY3Oo+3vKa+dT1y2dzwVwGlgMROIMr5wFvs6uuDKLiBQGvgV+A/5wb/bZ65Us1eeaTS65XsA5oADQBPgb8Lx7e6ZfM+tBee40EJJqW37glBdiyTSqug9ok2LTahEZh3OPbaJ3ospUp93fU147n79uAKraP+VrEekPHBSRCqq6x0thXRURqY7zy3sj0IU/r5NPX6/Un0tVLwE+f70A3J/lAvC7iEwFGrvfyvRrZj0oz23E+asnpRr8tVvrc0Sktoi8nmpzEM5fST5PVY8C+/jrtfP56wYgIm+ISM0Um4Lc333i2onITTi9i2+AB1T1XG64Xml9Lvd2X79ebUTkj1Sbg4Esu2bWg/LcfMAlIv1wyinvx7l344td9JSOAQNEJBb4EIgA+gC9vBpV5voEeNU99h8KDATGeDekTFEPaCwind2vxwA/qOpBL8bkERGpCnwPvKyq41K97bPX6wqfy2evl9tqoJy75zcGaIZTXXofToLK9GtmPSgPqeoFoCNOYjoCvAzc60P/caXJfYP2bpyKmxM4JbJvquqXXg0sc70CrMepulyB8xknezWizPE0zl+v24AYnGGXx70Z0FXoCRQEhorIqRRf7+Db1yujz+XL1wtVPQ7cDnTC+R04Ffi7qi4gi66ZrahrjDEmR7IelDHGmBzJEpQxxpgcyRKUMcaYHMkSlDHGmBzJEpQxxpgcyRKUMcaYHMke1DV5lohMw5mFOT2v48xEPR8oqKrZMt2OewLfxcATqrrlGo5PAu5S1e8zIZaSwC2q+nkmtDUQKKCqqWcuMSZN1oMyednzOLMwl8FZTgWgaYptI4Al7p9Pp3F8VukDrLmW5ORWBvg5k2IZjvNgZmYYBzzhnqfOmCuyHpTJs9xPxh8HEJES7s0H3RPoppT6dZYRkXzAS0D7a20jjfivhyuzGlLV8yIyHRgCPJlZ7ZrcyxKUMRlwr8F0eYjPPXz2KE4SEZylSh4DBuFMW3MCeElVP3EfXxAYibOMSRIwD3g+gyU0HgGOqep69/HhwE6c6ahGAeWBX3CWQxkB3AXEAz1Uda77mMtDfCISDSwA6gMdgD04S5B/4N43BhihquNTna+uO+auyW2qqktEAoE3cRJMCLAM6KOq6t6vtfvz1sWZ1udT979HovvzfQX8ISKDfH2aMJP1bIjPmKs3DOgLNAcqAitxElMTnF/AU0QkeR2qqTiJ7FacZU2ScJYAT++PwzuAH9PY/ibQGbgFZzhyLc7wYyNgFc4ik+l5EWfILwInuU0SEU+WHR8BzMCZ/LSMe9sb7hgfxpksVIEFIlLYfe/sG/f+NYEngGdI0VtyJ96DOGsiGZMhS1DGXL0JqjpfVVfj/DI+BQxx9yJG4fQsKotIFZweUWdVXeH+5fw4zrLgt6XTdmOcCTdTe8vdxiKcFVs3qupYVd0MTAAquHtraYlW1Qnu+F7CGTmpd6UP6S4KOQucV9V9IhIC9AOeVdUFqrpZVfvgDJM+DhQGigHxqhqjqr/g9Np+SdX0Rv5cQ8iYdNkQnzFXb1uKn88AMaqaPOty8to+wUAl988q8pelxPLj9KrSqrIrBRzy4Jwp90l5zpNpHHu52EJVT7hjCUxjvyup6j7HXPcwYrJ8gKjqEREZCkwVkdeA/wL/VtUVqdo5DIRdw/lNHmMJypirdzHV60vp7Bfg3jcCZ2gvpSPpHHOJtAsTPD1nWi6ksS35HKnjyuh3QvJ7twAHUr13AkBVh7jL9+/BWZrhZxH5p6q+nWJff3xkkT7jXTbEZ0zW2YTTUymgqttUdRtOQcO7QHql1vuAktkUHzjJq3CK11VSvZ8ygW0DEoCwFJ9nO/Aa0FREKorIJGCXqr6rqu2Ad3CKSFIqQTZWRhrfZT0oY7KIqqqIfAt8LCI9cYoD3sIprticzmF/4FTcZZcVQHcR+S/OEuT/x1+T0imgjohUUtVdIjIRGCMi54GtOPek7gVexfl8DwCIyEigEE5vK/UQXz2cVamNyZD1oIzJWl1xStG/wflFXRhnZoZj6ez/A061X3Z5Gece1VLgY5zeUMrhw2lAaWCTu/LvBeBL4COcSsKGQEfV/2/vjm0QBoIgik7sAmhhAyTqIKACUofQCzXQEHJAMgFtQOpg7YTEyUm3sv6rYLPR2XO3/tj+SbooK+aTshzxVl48liRFxFG5cbbVRWLsGBt1gUIiYlCuAz/bfnUep7mlRHGwPfaeBfVxggIKsf1V/qO69Z6ltaWmflU+nwRsIqCAeh6STvHXTd+Bu6Tn+uoEsIVPfACAkjhBAQBKIqAAACURUACAkggoAEBJBBQAoKQZ+TM/naKLQvEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(results.T, label='coffee')\n",
    "decorate(xlabel='Time (minutes)',\n",
    "         ylabel='Temperature (C)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's the final temperature:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "72.2996253904031"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coffee.T_final = get_last_value(results.T)\n",
    "T_final = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Encapsulation\n",
    "\n",
    "Before we go on, let's define a function to initialize `System` objects with relevant parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(T_init, r, volume, t_end):\n",
    "    \"\"\"Makes a System object with the given parameters.\n",
    "\n",
    "    T_init: initial temperature in degC\n",
    "    r: heat transfer rate, in 1/min\n",
    "    volume: volume of liquid in mL\n",
    "    t_end: end time of simulation\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    init = State(T=T_init)\n",
    "                   \n",
    "    return System(init=init,\n",
    "                  r=r, \n",
    "                  volume=volume,\n",
    "                  temp=T_init,\n",
    "                  t_0=0, \n",
    "                  t_end=t_end, \n",
    "                  dt=1,\n",
    "                  T_env=22)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we use it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "72.2996253904031"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coffee = make_system(T_init=90, r=0.01, volume=300, t_end=30)\n",
    "results = run_simulation(coffee, update_func)\n",
    "T_final = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "**Exercise:**  Simulate the temperature of 50 mL of milk with a starting temperature of 5 degC, in a vessel with the same insulation, for 15 minutes, and plot the results.\n",
    "\n",
    "By trial and error, find a value for `r` that makes the final temperature close to 20 C."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20.00135627897414"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "milk = make_system(T_init=5, t_end=15, r=0.133, volume=50)\n",
    "results = run_simulation(milk, update_func)\n",
    "T_final = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxU1fnH8U8SAiRA2JcQlrDIARSQVVDcV7TWrVoVUVvburVWLdWqLYrLT6u4VW1xqVYtLrW41X0pKhbZQQTCAwEBEwhhhxCyzvz+uBMcYhIGyMydJN/365VXMveee+cJZOaZc+65z0kIBoOIiIjEm0S/AxAREamKEpSIiMQlJSgREYlLSlAiIhKXGvkdQKw555oAw4H1QLnP4YiINHRJQDowx8yKw3c0uASFl5ym+x2EiIjs5Wjgy/ANDTFBrQeYMmUKnTp18jsWEZEGLS8vj7Fjx0LovTlcQ0xQ5QCdOnWiS5cufsciIiKeH1xy0SQJERGJS0pQIiISl5SgREQkLilBiYhIXIrpJAnn3MnAfcAhQD7wgJk96ZxrDDwO/ATvQtlDZnZvNedIAO4CfgU0Bp4Dfm9mZbUR444dO8jPz6e0tLQ2TlcvNGvWjC5dupCYqM8zIhI7MUtQzrmuwFTgMuAtYCjwoXNuNXAc4IBeQEvgA+dcrpm9UMWpfgWcCwwBioE3gFuBOw82xh07drBhwwYyMjJISUkhISHhYE9Z5wUCAXJzc9m0aRMdOnTwOxwRaUBi+ZE4E3jJzN4ws4CZzQE+A47CS1r3mNlWM1sNTAKurOY8lwGPmFmOmW0E7qih7X7Jz88nIyOD1NRUJaeQxMREOnbsyPbt2/0ORUTiSCAQZN2mArJzthGtZZti1oMys+mEVXBwzrXBu3P4RbwyF0vDmi8DBlRzqv5VtO3snGtjZlsOJsbS0lJSUlIO5hT1UnJyMmVltTKCKiJ1UHkgSG7+TlbmbmdlznZW5m5jVe52Cou894X7rh3NoT3b1vrz+nKjrnOuJfA2MAuYF9pcGNakEEit5vDmVbQl1P6gEhSgnlMV9G8i0nCUlgVYm7cjlIy2sTJ3O9+u20FJ6Q9Ll7Zu0YT+PdvSpUPzqMQS8wTlnOuDdw1qKTAWqOiyhHddUoGCak6xq4q21NBeRESqUFRSxur1O1gV1jNas34HZeU/HLLr0DqFXl1a0SujJb26tKJnRkvapDWNanyxnsV3DF5ymgzcamZBoMg5l4c3SSI31LQvew/jhVsaavu/sLbrzWxb1AKvw9atW8cZZ5zBF198wfbt2znxxBOZM2cOO3bs2PNzWlqa32GKSJQVFpV6iSisZ5SzYSeBKi4fZbRvRs+MimTUkp4ZrUhr1jjmMcdyFl8v4B3gNjN7rNLuF4HbnXOL8IbwxgOPVnOqF4HxzrlP8XpTd4S2SRU6d+7MggULADTRQaSBCAaDrNu0i6xvt7BszRaWrd7C2g07qTyXITEBundqsVfPqEfnNFKbJvsTeCWx7EFdC7QA7nXOhd/j9AQwAXgQWII3s/ApvF4WAM65AuBKM5sS2t4RmIE3vPda6PgGJScnh7PPPpvf/va3TJ48mdLSUn7zm9+QmJjI5MmTKSkp4eqrr+akk07a01OqTnl5OePHjyc/P5+nnnqKZs2axfA3EZGDVVRcxorvtrFszRayVm9h2eqt7Cws2atNo6QEMtPTvJ5Rl5b0ymhJ9/Q0mjaO35rhsZzFdyNwYw1Nrg19VXVs87CfA8Dtoa+omvjMTOZmbYj20+wxrF9Hbv/FyIjb79y5k8WLFzNt2jQ++eQTxo8fz1lnncWnn37K9OnT+c1vfsOIESNqPEcgEODWW29l48aNPP3006SmVjc3RUTiQTAYJH/rbrJWb8FWbyFrzRa+XbeDQKWxulYtmtAvsw19u7emb2YbenVpRZPkJJ+iPjDxmzolItdccw2NGzdm1KhRlJeXc+mll9K4cWOOP/54ysvL91n9YeLEicyfP5/3339fyUkkDpWUlrMyZ3tY72gLW3futfAsiQnQM6Mlfbu39pJSZhs6tqn793MqQdVgf3ozfmnVqhUASUneJ6MWLVoA7ElM+7qBbv369ezcuZM5c+Zw7LHHRjFSEYnE1h1FLA1dO8pavYWVOdspKw/s1aZFajKuexv6ZnoJ6ZCurUlpUv/ezuvfb9TAHOwnpCeeeIJ3332XCRMm8O6779K8eXTuZxCRqhUUlvDNys0syt7I1ys28d2GnXvtT0iAbp1a7DVcl9G+eZ3vHUVCCaqBS05O5pJLLuHtt9/m/vvv5847D7qkoYjUoKi4jKWrt7BoxUa+zt7Eqpxte031btI4iX7dvWG6fplt6NO9Nc1T4mNWXawpQQmJiYncddddnH/++Zx++umMHBn/Q5sidUVpWYDla7fuSUi2ZsteN8I2SkqgX/c2DOrdjoGHtKdPt9YkN9LKAaAEVWd16dIFM9vzOC0tba/HwJ7HFd/D21Ru369fPxYvXhztsEXqvfJAkFW521i0YhOLsjex5NvNFJd8XyYoIQF6d221JyH1z2xD03p4/ag26F9FROQgBINBvtuwk69XbGJR9ka+WbmZXbv3Xk+uW6cWDOzdjoG92zOgV1uap8a+KkNdpAQlIrKfduwqYd6yDczLymdR9sYfTPvu2CaVQYe0DyWldrSOcs26+koJSkRkHyp6SbOXbmDO0jyWrd6y18SG1i2aMLB3ewYd0o4BvdvRqa2qsdQGJSgRkSqUlgVYsmrTnqSUt/n7VX4aJSUwsFc7hvfryGDXgS4dGsa071hTgqokGAzqD62SaK2WKRJvthcUMzdrA3OWbmC+5bO7+PuFOtOaNWZYv46M6N+Jwa593BRUrc+UoMIkJyeze/dulfyppLS0lEaN9Kci9U8wGGRN3k7mLM1j9pI8bO3WvSp+Z6anMby/l5QO6daapER9eI0lveuE6dChA7m5uWRkZJCSkqKeFF4x2Q0bNtCyZUu/QxGpFaVl5XyTvZnZS/OYszSP/K279+xrlJTIwN7tGNG/I8P6d6JjG31Y9ZMSVJiKhfvWrVtHaWnpPlo3HM2aNaNdu3Z+hyFywLbtLGZuVh6zl25ggeVTFHZfUqsWTRjeryPD+3fk8D4d6mVNu7pK/xOVpKWlaYVZkXpgZ2EJMxatZ/rCHL7J3rTXrLuenVt6Q3eHdqJ3l1YkauguLilBiUi9UVhUyuwleXyxMJcFlr+npFCjpESG9GnvDd3160T71ik+RyqRUIISkTqtuLScuVkbmL4glzlL8ygp85amSExMYHCf9hwzOIORAzo32IKrdZkSlIjUOWXlARYu38gXC3KYuThvr+ngh/Zsy9GHZ3DUwM60atHExyjlYClBiUidUB4IsmTVJr5YkMuMRevYWfj9RKbeXVtxzOEZjB6UoeG7ekQJSkTiVjAYxNZu5YsFuXy5MHevmnddO7bg2MEZHH14Bp3ba6HN+kgJSkTiSjAY5Nt1O/hiQQ7TF+budZ9Sp7apHH14BscM7kJmumbb1ndKUCISFzZv380nc9by2bwccvIL9mxvk9Y0lJQyOKRrK91A34AoQYmIb8oDQRZYPh/OXM3spRsIhG5WapHamNGDOnPM4Az692ir+5QaKCUoEYm5zdt38/HstXw0aw0bQ0N4iYkJjBqQzilHdOfwPu1plKRlzxs6XxKUc24E8I6ZdQg9LqjUpBHQBMgws3VVHL8WaAtU3Buea2YuiiGLyEEqDwSZt2wDH361hrlZeXsqO3Rsk8qpI7tz4vButNHCfhImpgnKOZcAXAFMCt9uZs3D2jQCpgGfVZOc2gEZQJqZ7YpuxCJysPK3FvLJ7LV8PGsNm7YXAZCUmMBRA9M5bWR3BvZuryE8qVKse1ATgTOAu4E/VtPmZiAZuL2a/UOBFUpOIvGrvDzAnKwNfDhzDfOXbdjTW0pv14xTj+jOCcO70rqFektSs1gnqMlmNsE5d1xVO51znYFbgaPMLFDNOYYAic652UAPYD5wvZllRSNgEYnchi2FfDxrDR/PXsuWHV5vqVFSAqMHdOaUkd0Z0KudeksSsZgmqKqG7Cq5AfjAzBbW0KYcmA38AdgCTADed871N7PCGo4TkSgoKw8we0keH85awwLL37PgX0b7Zpw6MpMThnWlZXOVHJL9Fzez+JxzScBlwNia2pnZ/ZWOuwW4Bm/ob3rUAhSRveRt3sVHod7StlCFh+RGiRw1sDOnjuzOoT3b6p4lOShxk6CAI0PfP62pkXPuemCumX0Z2pSE93sURTE2EQlZ+u1mXp+WzawleXu2de3YnFNHZnL80K6kNWvsY3RSn8RTghoJzKzh2lOFTGCcc+5HwDbgz8AKvGtRIhIFgUCQWUvyeH3aCpat2Qp4vaWjD8/g1JHd6ZfZRr0lqXXxlKAygSqvUYXuk7rSzKbgXXt6EFgANAM+B840s/KqjhWRA1dSWs60ed/xxmfZ5G70Js42T0nmjKN6cMboHpqJJ1HlS4Iys8+AVpW2XVtD++ZhPxcB14a+RCQKCgpLeG/Gav7z5ao915c6tE7hrGN7cfKI7qQ0iafPtlJf6a9MRPbI31LIW9NX8tHMNRSVeIMSPTNacu5xvRk9qDNJKj8kMaQEJSKsyt3O69Oymf517p6CrYP7tOfc43sz6JD2ur4kvlCCEmmggsEgX6/YyNRp2SxcvhHwCrYeN6QL5xzXm54ZLX2OUBo6JSiRBqa8PMCXX6/j9WnZrFq3HYCmjZM4ZWR3zjq6Fx3apPocoYhHCUqkgdhdXMbHs9fw1ucr96xS26pFE84c3ZPTj8ykearuX5L4ogQlUs9t21nMO1+u4t3/fUvB7lLAK0N0znG9OX5oVxonJ/kcoUjVlKBE6qmCwhKmTsvm7emrKCn1ZuT1y2zDOcf15ohDO6loq8Q9JSiReqaopIz/TF/F1GnZ7Ar1mI44tBPnHt+b/j3a+hydSOSUoETqidKyAB/NWsOrHxtbQzfXDjqkHZee3p8+3Vr7HJ3I/lOCEqnjAoEgXyzIYcqHy8jb7K0407trKy47vR+H9+ngc3QiB04JSqSOCgaDzMnawIvvZbF6/Q4AunRozrgx/Rg1IF0310qdpwQlUgctWbWZ599dStbqLQC0a5XC2FMdxw/tqnJEUm9ElKCcc2OAMcAwoAPeqrZ5wBzgnVDxVxGJslW523nhvaXMW5YPQFqzxlxwUh/GjMrUdHGpd2pMUM65i4A7gDbAJ8D7wGa8RQLbAYOAfznnNgL3mNlLUY1WpIFat6mAKe8v44uFuQCkNGnEOcf24qxje5HaNNnn6ESio9oE5Zz7ENgFXAV8Xt1Cgs65RLze1dXOuZ+Z2clRiVSkAdq8fTevfLycj2etoTwQpFFSImcc1YPzTzyEls2b+B2eSFTV1IOaYGaz9nWCUOJ6F3jXOTey1iITacB2FpYw9b8r+M+X31JSWk5iApw8ohsXnuLo0Fq18qRhqDZBRZKcqjhm5sGFI9KwFRWX8fb0Vbw+bQW7isoAOHJgOpec1o+uHVv4HJ1IbO3rGlQvYCJwm5mtCdv+D6ApcHP4dhE5MIFAkI9nr+GfHyzbs4Lt4Ye0Z9zp/XSTrTRYNV2DOgT4H7ABaFZp9zTgRmCWc26UmX0bvRBF6reVOdv429RF2NqtAPTp1opLT+/PoEPa+xyZiL9q6kHdBXwFnGdmZeE7zOx559zLeNee7gTGRS9EkfqpsKiUKR8s450vVxEIQpu0Jlzx48M4+vAM3WQrQs0J6hjgzMrJqYKZlTjn7gRejkpkIvVUMBhk+sJc/v72YrbsKCYxAX58TE/GntpXU8ZFwtSUoFoAW/dxfA6gdaFFIpSTv5PJry/i6xWbAOjbvTVXnzdIy6uLVKGmBJUFjARW1dBmFLC6NgMSqY+KS8t57ZPlTJ2WTVl5gBapyVz+o0M5aXg3rcskUo2aEtQzwP8552aY2erKO51zPYH7gMeiFJtIvTBnaR5PvvENG7Z4lcZPHtGNy87orxttRfahpvugnnLOHQcsdc79G6/u3nagNTAcOBev9NFDMYhTpM7J31rI029+w8zFeQBkpqdxzXmD6Nejjc+RidQNNd4HZWYXO+cuBn4O3A60wqvFNwu4xMxeP5Andc6NwCsy2yH0uAmwEygJazbDzE6p4tgEvBmGvwIaA88Bv69uModIrJWWBXj7i5W8/LFRXFJOSpMkLj61H2eO7qFK4yL7YZ/VzEMFYGulCGwouVwBTKq0awCwxcw6RXCaX+H13oYAxcAbwK14091FfPXNyk38beoivtuwE4DRgzrzi7MOo23LFJ8jE6l7arpR9zbgITPbHcmJnHPNgN+ZWU2JYiJwBnA38Mew7UOBhZE8D3AZ8IiZ5YSe9w7geZSgxEfbdhbz3DtL+O/c7wBIb9eMq84dyBCnFW1FDlRNPahiYLFz7lXgdTObW1Uj59xg4BLgfODxfTzfZDObELq2FW4I0ME5twjoCHwBXG9muVWcoz+wNOzxMqCzc66NmW3Zx/OL1KryQJAPZ67mhfey2LW7lORGiZx/Yh/OO7631mcSOUg1TZKY5Jx7DbgZ+Mw5txsvMWwCEoD2wGF414FeAI6parZfpXOuq2bXLryySncCpcBf8IbuRlTRtjlQGPa44udUQAlKYib7u238derXrPhuGwBD+nbgqnMGkt6ucmUwETkQ+5oksQa4xjl3M3AC3lBcRyAAfA38Gfg00mHAGp7nxvDHzrkbgY3Oua5m9l2l5ruA8AH9irUHCg4mBpFIFRaV8sJ7Wbw341uCQWjXsim/PHsAowakq0SRSC2KaMl3M9sJvBX6qnUVJZPMLCu0qXHoe1EVzZcCDq/HBdAXWG9m26IRm0i4rG+38NDL88jbXEhiYgJnH9uLi05xpDSJ6KUkIvshXl5VA4FhoSntAI8C75rZxiravgiMd859itebuiO0TSRqysoDvPKx8donywkEoWfnltxw8RAy09P8Dk2k3oqXBHUF3nWnbLyY3sWbTg6Ac64AuNLMpgCT8YYZZ+AN770GTIh1wNJwrNtYwIMvzWP52m0kJMB5x/dm7Gn9SG6ke5pEosmXBGVmn+Hd9FvxeDMwtob2zcN+DuDdNHx7FEMUIRgM8tGsNTz91mKKS8pp3zqFGy4awoBe7fwOTaRBiJcelEhc2V5QzGP/WsisJV6ZomMHd+Gq8wbSPEXLYYjESsQJyjk3DG8V3T7A2cBPgVVm9kaUYhPxxdysDTz66gK27SymWdNGXHXeII4b0sXvsEQanIgG0Z1zpwKfA7uBQ/Fm2aUCrzrnLo9adCIxVFRSxt+mfs3EZ2aybWcxh/Vqy1/GH6/kJOKTSK/y3o1XxugKoAzAzO4CrsO7kVekTsvO2cYND3/OezNW0ygpgcvP6M/dVx1Fh9ap+z5YRKIi0iG+/sBHVWz/CHi49sIRia3yQJDXp61gygfLKA8E6dqxOb+7eCi9urTa98EiElWRJqhcYDA/XF33BGBtrUYkEiMbthTy8MvzWbJqMwA/Gt2Dy390KE1UQ08kLkSaoO4FngqtopsEnOacywSuBW6IUmwiUREMBvlsfg6TX19EYVEZrVs04bcXDmZo345+hyYiYSItdfSccy4fuAWv/NB9eFXELzOzf0cxPpFaVVBYwt+mLuKLhV6h/JGHdeLX5x+u5ddF4lBECco5dxMwxcxGRzkekahZlL2Rh1+az6btRTRtnMQvzx7AySO6qcCrSJyKdIjvNmBqNAMRiZbSsnJefH8Zb36eTTAIrltrbhw7hM7tmu/7YBHxTaQJ6l3gN865iWa2NZoBidSmNXk7mPTPeaxev4PExAQuOrkPF5zUh6Qk1dETiXeRJqjewIV4SWoH3g27e5hZ59oOTORgTV+QyyOvLqCktJz0ts248eIh9M1s43dYIhKhSBPUk6EvkbhXHggy5YMsXvt0BQDHD+3CVecOJLWp6uiJ1CWRzuL7e7QDEakNu3aXMmnKPOZmbSAxMYErfnwoZ47uqYkQInVQpLP4Xqppv5ldXNN+kVjIyd/J3c/OJndjAS1Sk7l53HAG9Wnvd1gicoAiHeIrr+K4nsAAvIUGRXw1N2sDk/45l11FZXTv1II//vwIOrVt5ndYInIQIh3iG1fVdufcrUCvWo1IZD8Eg0GmTsvmhfeWEgzCqAHp3HDREFKaaKkzkbruYF/FLwGL8JZsF4mpopIyHnt14Z6qEGNP68sFJ/YhMVHXm0Tqg4NNUJcDO2shDpH9kr+1kHuem82q3O2kNEnihouGMmpAut9hiUgtinSSxHogWGlz89DX72o7KJGaLFm1mXufn832ghLS2zbjtp+PoHunNL/DEpFaFmkP6o+VHgeBEmCOmVnthiRSvfdnfMuTb3xDeSDI4X3ac9O4YbRIbex3WCISBZEmqI7Aw2a2VwUJ51yac+5+M7up9kMT+V5pWYCn3vyGD75aDcDZx/bi8jP6q2SRSD1WbYJyzqUDLUIP7wKmOec2V2o2GPgNoAQlUbN1ZxH3PT+Hpd9uIblRIr8+/3BOGNbV77BEJMpq6kEdCbzG99ee/he2LwhUTJX6R+2HJeLJztnGPc/NZtO23bRJa8ptPxtBn26t/Q5LRGKg2gRlZlOdc72BRGA5MArYFNYkCBSYWf7+PqlzbgTwjpl1CD3uADwKnIiX+N4Hfltd5XTn3FqgLd8nz1wzc/sbh8S3z+fn8JdXF1BSFqBv99bccvkI2qQ19TssEYmRGq9BmdkqAOdcsplVriZB2L7SSJ7MOZeAd8/UpEq7ngG2Az2AZOBF4AngByWUnHPtgAwgzcx2RfK8UreUB4K8+N5Spk7LBuDkEd24+ryBJDdK8jkyEYmlSCdJtHHO/QHoD1S8SyQATYDDgEjXMJgInAHcTWhmoHMuEQgAEysSjnPuaeDxas4xFFih5FQ/Fewu5YF/zmX+snwSExP45VmHccZRPVTsVaQBijRBPQUMA/4D/DL0uDcwGvj9fjzfZDOb4Jw7rmKDmQWAsyu1OxtYUM05hgCJzrnZeD2u+cD1Zpa1H3FIHPpuw07ufnYW6zbtokVqY265bDgDerfzOywR8Umkc3RPAC4zs2uAJcCzZnYqcD/eZIqImNm6fbVxzo3HS1A3V9OkHJgNnAt0x0tk7zvnUiONQ+LP/GX5jP/LF6zbtIvM9DQevuFYJSeRBi7SHlQTYEXo5yy8YbZ5wLPAF7URiHMuGXgMOBM4wcyWVdXOzO6vdNwtwDWhmKbXRiwSW9MX5PLgS/MoDwQ5alBnrv/pYJqq2KtIgxdpD2o5cFTo5yzgiNDPqaGvg+KcawF8DAwHRpjZwhraXu+cGx22KQkv0RYdbBwSex/OXM0DU+ZSHghyznG9uXncMCUnEQEi70FNAv7hnEsCXgW+ds4FgZHUTq/lFbxkebSZFe6jbSYwzjn3I2Ab8Ge83t38WohDYuj1adk8984SAMaN6cf5Jx6iyRAiskdEPSgzewE4CcgK1d47C2/m3lcc5FIbzrmBwOnACCDfOVcQ+soJa1PgnBsbevgHYCbetad8vIUTz6xuGrzEn2AwyAvvLd2TnK46ZwAXnNRHyUlE9pIQDFYuUv5Dzrm3gJuruy5UlzjnMoFvP/30U7p06eJ3OA1OIBDkyTcW8d6M1SQmJnD9hYM5fqjKFok0VDk5OZx44okAPcxsdfi+SIf4RuNVLxc5YGXlAR59dQGfzcuhUVIiN186jJGHaQ0nEalapAnqIeB559xDwCpgr6rmZra8tgOT+qWktJz7X5zLrCV5NG2cxB9/fgSDDmnvd1giEsciTVB3hb5XzOSrGBdMCP2sGjRSrcKiUu55bjaLsjfRPCWZO345Etc90uIjItJQRZqgDolqFFJv7dhVwsRnvmL52m20SWvCnb86ku7pWv1WRPYtogRlZisBnHN9AAd8CrQzs7VRjE3quM3bdzPhqa9Ym7eTjm1SuevKI0lv18zvsESkjogoQYVupP0nXpWHANAHeDA0I+6MSEoYScOSt3kXf3pyBnmbC+nasQV3XTmKti1T/A5LROqQSCtJTAJa4t0kWzFB4kagAHik9sOSumxN3g5ufnw6eZsL6d21FfddO1rJSUT2W6QJ6kfA78KH9MzsW+A6vEUGRQBYvnYrtzzxJVt2FHNYr7bcc9WRpDVr7HdYIlIHRTpJogVQXQkiFU4TAL7J3sRdz85kd3E5w/t35OZLh9MkWRM8ReTARNqD+hCYEKrFBxB0zrUFHsAr8ioN3Owledz+9FfsLi7n2MFduPXyEUpOInJQIu39XAe8AWzCq17+AdAFr0jr5VGJTOqMz+bn8PDL8wkEgow5MpOrzhlIYqLq6onIwYl0mvl6YKRz7mS8Zd8b4S278UFoRVxpoN6b8S2TX19EMAjnn3gI48b0U9FXEakV+3v9aBVeD6ocWKzk1LC99ulyXngvC4DLzujPT07Q/dwiUnsivQ+qFfAvvBl75XgljhKdc28Al0awhpPUI8FgkOffXcrUadkkJMDV5w1izKhMv8MSkXom0kkSTwLt8NZsaoq3BPwooBe6D6pBCQSC/HXqIqZOyyYpMYHxY4cqOYlIVEQ6xDcGOM7Mwletne2cuxJvwsSvaj0yiTvBoLeW0wdfraZxo0T+cNlwhvfv5HdYIlJPRZqgNgGtqzl+V+2FI/FsyofLeG/GapIbJTLhFyO1XIaIRFWkCeoW4Cnn3C3Al0AZMAR4GHgsVEQW0NpQ9dXbX6zk1Y+Xk5iYwE3jhik5iUjURZqgXg59f4W914ICuA+4F60NVW/9d+5ann5rMQDXXXC4VsEVkZjQelBSo1mL1/PoqwsBuOLHh3Hi8G4+RyQiDcV+rQclDcs3Kzfx5xfnEggEueCkPpx9bC+/QxKRBiTS+6AGAg8Bh+JNMd+LmWn97nomO2cbd/19FqVlAcaMyuSS0/r6HZKINDCRDvE9jzdb749AUfTCkXiQu7GAO57+it3FZRx9eAZXnjtQ5YtEJOYiTVB9gOFmtjSawYj/Nm7dzZ+enMH2ghKG9O3ADRcNIUmFX0XEB5FWkvgf3vCe1GPbC4qZ8NQMNm7dTb/MNtxy6XCSG0X6JyIiUrsi7UH9Evifc24MXsHYvYrEmv2NDgAAABYASURBVNn/7c+TOudGAO+YWYfQ48bA48BP8Gr9PWRm91ZzbAJwF171isbAc8Dvzaxsf2KQvRUWlXLHMzPJyS8gMz2NCVccQdMmWotSRPwT6TvQn4B04ChgcKV9QSCiBBVKLlcAkyrtmgg4vNp+LYEPnHO5ZvZCFaf5FXAu3o3CxXjrVN0K3BnRbyI/UFJazj3PzSb7u210apvKxF+NonmqlmkXEX9FmqAuBM41s7cO8vkmAmcAd+NNuKhwGXC5mW0FtjrnJgFXAlUlqMuAR8wsB8A5dwfeJA4lqANQXh7ggX/OZVH2Jlq3aMJdVx5Jm7SmfoclIhLxNaitQG2UMJpsZkOBuRUbQkt5pAPhEzCWAQOqOUf/Ktp2ds5pqvt+CgaDPP7a18xcnEezlGTuvPJIOrVt5ndYIiJA5D2om4G/OOd+B6wESsN3mllJJCcxs3VVbG4e+h6+plQh3sKIVWleRVtC7bdEEod4yenZ/yzhkzlraZycxO1XjCQzPc3vsERE9og0QT0ItAUWVLP/YOrvVVRDTwnblgoU1NC+cltqaC9V+Pd/V/Dm5ytplJTArZcPp18PdUBFJL5EmqAuiVYAZrbVOZeHN0kiN7S5L3sP44VbGmr7v7C2681sW7RirG/e/2o1L7yXRUIC3HjRUIb27eh3SCIiPxBpLb5PK352zrUACswsWMMh++tF4Hbn3CK8IbzxwKM1tB3vnPsUrzd1R2ibRGD6wlz+NvVrAK4+dyBHD87wOSIRkapFfBemc+5G59x6vAkTPZxzzznnHnHO1cbNMhOAxcASYA4wFZgc9twFzrmxoYeTgdeAGcAKvB7VhFqIod6bvyyfh16aRzAIl4zpy5gje/gdkohItSItFvs74DrgJr5PHG8BfwNKQtsjZmafAa3CHhcB14a+qmrfPOznAHB76EsitGz1Fv7v+dmUlQc565heXHBin30fJCLio0h7UL8ErjKzFwlVkTCzN/HuSbo4SrFJLVmzfgcTn5lJcUk5Jwzrys/PPFTFX0Uk7kWaoLrh3W9U2WpA07/iWN7mXUx4agYFu0s54tBOXHfB4SSq+KuI1AGRJqgFwAVhjysmSFxF9VPPxWdbdxQx4cmv2LKjmAG92nHTuGEkJan4q4jUDZFOcBiPVx/vGLwFCyc65/riVXUYE63g5MAVFZdxx9MzWb95F727tOSPPx9B4+SDuV1NRCS2Ivo4bWZf4a0JtQB4F2gNfAb0M7PpUYtODkgwGOQv/1rIqnXbSW/XjDt+OYrUpsl+hyUisl+q7UE55yYAk8ysEMDMNrB3gVeJU298ls30hbmkNEnijz8bQcvmTfwOSURkv9XUg7qd7+vkSR0x3/J5/l2vCMcNFw2lWyfV1xORuqmmBKWpXnXM+k27eODFuQSCcOHJjlED0v0OSUTkgO1rkkQX59w+Fwcys7W1FI8coN3FZdzz3CwKdpcyon8nLjrF+R2SiMhB2VeCmrOP/Ql4U841PcxHwWCQR19ZwJq8nWS0b86NFw/RvU4iUuftK0EdD2yORSBy4P793xX8b9E6Upo04rafjaBZimbsiUjdV1OCCgLLzCw/VsHI/pubtYEX388CYPzYoXTt2MLniEREaocmSdRh6zYWMOmfcwkG4eJT+zLi0E5+hyQiUmtqSlDPA7tjFYjsn8KiUu5+bja7isoYeVgnfnqSqpOLSP1S7RCfmf0sloFI5AKBII+8soDvNuyka8fm3HCRJkWISP2jyqF10GufLuerb9bTrGkjbvvZESpjJCL1khJUHTN7aR5TPlxGQgKMv2QYGe1V7ENE6iclqDokJ38nD04JLdl+Wj+G9evod0giIlGjBFVHFBaVcs9zsyksKuPIgemcf+IhfockIhJVSlB1QCAQ5KGX5pOTX0D3Ti24/sIhWrJdROo9Jag64JWPjVlL8miWksxtPzuClCaRrjMpIlJ3KUHFuZmL1/PyR0ZiAtx0yTDS2zXzOyQRkZhQgopj323YyUMvzQNg3On9GdK3g88RiYjEjhJUnCrYXcrdz85id3E5owd15rzje/sdkohITMXFxQzn3FjgyUqbU4BPzeyUSm1PAD5m7zJMfzazu6IbZewEAkEenDKPdZt2kZmexm9/OliTIkSkwYmLBGVmU4ApFY+dc4OBj4DfV9F8CPCamV0Yo/Bi7qUPlzE3awMtUpO57WcjaKpJESLSAMXdEJ9zLhkvWd1hZl9X0WQosDC2UcXOjEXrePWT5d6kiHHD6NRWkyJEpGGKx4/m1+IN3/21mv1DgPbOuavxlgR5FfijmRXHKL6oWbN+Bw+/PB+Ay390KIf30aQIEWm44qoH5ZxrjDesd4eZBavY3wjIAd4A+gEnACcBdf76U0FhCfc8N5uiknKOHdyFs4/t5XdIIiK+irce1GlAAHi3qp1mVgacGLYp2zl3D/Bn4Kbohxcd5YEgD0yZx/rNu+jZuSW/vmCQJkWISIMXVz0o4CzgX2YWqGqncy7DOTcp1NOq0Bgoikl0UfLqx8b8ZfmkNWvsTYpoHG+fG0REYi/e3glHAn+qYf9mYCxQ6Jy7E+gB/BF4NgaxRcWyNVt49ZPlJIQqRXRok+p3SCIicSHeelCZwLrwDc65sc65AgAzKwLGAMfgJasvgNeAh2IbZu0oKi7joZfmEwgEOfvY3gzq097vkERE4kZc9aDM7AdzqivfI2VmC4HjYhhW1Dz7nyWsD92MO25MX7/DERGJK/HWg2ow5mZt4P2vVtMoKZEbLx5CcqMkv0MSEYkrSlA+2F5QzKOvLgBg3Jh+9Ojc0ueIRETijxJUjAWDQR5/bSHbdhZzWK+2nKX7nUREqqQEFWOfzvmOmYvzSGnSiBsuHEJSou53EhGpihJUDOVt3sVTb34DwFXnDtCUchGRGihBxUh5IMgjryxgd3EZRw5M5/ihXf0OSUQkrilBxcibn2WzZNVm2qQ14dqfHK5SRiIi+6AEFQOrcrfzzw+yALjup4NJa9Z4H0eIiIgSVJSVlJbz4EvzKCsPcvqRmQzt29HvkERE6gQlqCh78f0s1ubtJKN9M372o0P9DkdEpM5Qgoqir1ds5M3PV5KYmMCNFw/V0u0iIvtBCSpKCnaX8sgrXrWIC0/qQ59urX2OSESkblGCipIn31jEpm276dOtFeef1MfvcERE6hwlqCiYvjCXz+bl0KRxEjdePJRGSfpnFhHZX3rnrGWbt+/mr//+GoArzjyUjPbNfY5IRKRuUoKqRYFQtYiC3aUM7duB00Zl+h2SiEidpQRVi96b8S0Ll2+kRWpjrvvpYFWLEBE5CEpQteS7DTt57j9LAPj1+YNok9bU54hEROo2JahaUFYe4KGX5lFSFuCEYV05cmBnv0MSEanzlKBqwSsfGdk52+nQOoUrzxngdzgiIvWCEtRBWrZ6C699upyEBLjhoiGkNk32OyQRkXpBCeog7C4u46GX5hMIwrnH9eawXu38DklEpN5QgjoIf397Mes37yIzPY2xp/X1OxwRkXpFCeoAzV6ax4cz19AoKZHfjR1KcqMkv0MSEalX4qa8tnPu58CTQHHY5mvN7PlK7boBfwdGAvnAb8zsvZgFCmwvKOaxVxcCcOnp/chMT4vl04uINAhxk6CAIcCDZvaHfbR7BfgKOAMYDbzpnDvczFZFO0CAYDDIY/9ayLaCYgb2bsdZx/SKxdOKiDQ48TTENxRYWFMD51wfYBgwwcxKzOy/wNvAFTGID4BPZq9l1pI8Ups24rcXDiYxUdUiRESiIS56UM65JGAgMM459xBQCDwD/NnMgmFN+wNrzWxX2LZlwIhYxFlYVMrTb30DwFXnDqRD69RYPK2ISIMUFwkKaA/MBZ4HzgX6AW8BO4C/hrVrjpe8whUCMcsU7VqlMqBXW44b0iVWTyki0iDFRYIyszzg2LBNC51zjwHnsXeC2gWkVDo8FSiIboShJ2qazF9vOiEWTyUi0uDFxTUo59yhzrmJlTY3BooqbVsKdHPOhSepvqHtIiJSj8RFDwrYBvzOOZeDN4V8MHAd8OvwRmZmzrmvgXucc7cARwJnAaNiHK+IiERZXPSgzCwX+DFwJd51p6nAXWb2b+fcWOdc+BDeeXjXqPLxJlJcYWaLYx2ziIhEV7z0oAhNGR9WxfYpwJSwx98BY2IYmoiI+CAuelAiIiKVKUGJiEhcUoISEZG4FDfXoGIoCSAvL8/vOEREGryw9+IfLAnREBNUOsDYsWP9jkNERL6XDqwM39AQE9Qc4GhgPVDucywiIg1dEl5ymlN5R0IwGPxhcxEREZ9pkoSIiMQlJSgREYlLSlAiIhKXlKBERCQuKUGJiEhcUoISEZG4pAQlIiJxSQlKRETiUkOsJHHAnHODgMnAQGAV8HMz+8Hdz/HAOXcycB9wCN7ijg+Y2ZP+RlUz51wrYBEwwcz+4XM4VXLOpQN/A44HioCnzOxP/kZVPefcSOAvgAM2AveZ2TP+RrU359wI4B0z6xB63Bh4HPgJXrWXh8zsXh9D3EsV8XYAHgVOBBKA94HfmtlW/6L0VI610r7HgAFmdlzMA4uQelARCr1o3gJeBVoB9wAfOefSfA2sCs65rnirEt+NF+tFwL3OuVN9DWzfJgMZfgexD2/hlcnqCIwELnPOXexvSFVzziXixfsXM2uJ93fweOiDlu+ccwnOuV8AHwGNw3ZNxEuovYDheP/Gl/oQ4l5qiPcZoAzogfeBsDXwROwj/F4NsVbsHwNcHfPA9pMSVOSOA5LN7BEzKzWzV4AlwE/9DatKmcBLZvaGmQVCvbzPgKN8jaoGzrnLgDTgG79jqY5z7gigJ3CdmRWZ2bd4fxfTfA2seq2BDkCCcy4BCOK9kZb4GtX3JuK9Sd5daftlwD1mttXMVgOTgCtjHFtVfhBv6ENAAJhoZrvMbBvwNDDanxD3qO7fFudce+AxvJGAuKYEFbn+QFalbcuAAT7EUiMzm25mV1U8ds61wSuQu8C/qKrnnOsB3A783O9Y9mEoXgK9wzmX65xbCZxjZut9jqtKZrYZb6jseaAUrxjnrWZW+e/YL5PNbCgwt2JDaJg3HVga1i5eXmc/iDf0AfBsM8sOa3c2/r/WfhBrmGfxhv9XVrEvrihBRa45UFhpWyGQ6kMsEXPOtQTeBmbhDffEFedcEvBPYLyZxfsiXRWJvhSvJ3UuMD7Oh/iKgIuBFLze3u3OuVP8jKuCma2rYnPz0Pfw11pcvM6qiXcvzrnxeAnq5uhHVL3qYnXOXQME4u06ZHU0SSJyu/Be5OFSgQIfYomIc64PXlJaCow1s4DPIVXlT4CZ2et+BxKBYmCHmd0Revy1c+4ZvET1km9RVe9c4Cgz+33o8efOub/jDZd95F9YNdoV+h7+Wovr1xmAcy4Zb9jsTOAEM1vmc0g/4JzrB/weGOF3LJFSDypyS/Eu3Ibry95DEXHDOXcMXq/pTeAnZlbkc0jVuRD4iXNum3NuG95Qzl+dc3/1Oa6qLANSQxNmKsTzh7yuQJNK28rweoBxKTTzLY+9X2tx+zoDcM61AD7Gm9AxwswW+hxSdc7Fuya5IvRauxcYHfo5LsXziyveTMO72HwD3rj+eXjTzd/wNaoqOOd6Ae8At5nZY37HUxMz6xv+2Dm3EHgkTqeZf4w3VftB59zv8N5EryB+Z0N9hDd781d4F+6HAL8EfuFrVPv2It5Q5CK8Ib/xeNO449UreB/2jzazypcB4oaZ3YM3+xgA59z1wNmaZl4PmFkJMAYvMW0BbsP7z93oa2BVuxZogffmVBD29We/A6vLQr3QY/GuP60HPgDuN7OpvgZWDTNbgvep+UpgG94w5B/MLO6uRVYyAViMN0t2Dt4tE5N9jagazrmBwOl4w2b5Ya+1HJ9Dqxe0oq6IiMQl9aBERCQuKUGJiEhcUoISEZG4pAQlIiJxSQlKRETikhKUiIjEJd2oKw2Wc+4feJWzqzMRrwr8NKCFmcWk3E6oPuH/gEvNbPkBHB8EzjSzd2ohlvbAyWZ20KWcQnXqmpnZxIM9lzQM6kFJQ/ZbvMrZ6XiFVMG74bJi2yRgRujnXVUcHy3XAV8fSHIKScerelEb7se72bc2PAZcGqoRKbJP6kFJg2Vm24HtAM65dqHNG6uoqh6zKuvOuabALcAJB3qOWq4Kn1BbJzKzYufc88CtwOW1dV6pv5SgRGrgnDuOsCG+0PDZRXhJxOGtt3MJXpXoccAO4BYzezF0fAvgQbzly4PAf/GWA69u6YYLgW1mtjh0fCbwLfBj4CGgC/AJXv2/SXjVs9cD15jZR6Fj9gzxOec+Az4HBgGnAN8BD1Qst+CcWw1MMrPHKz3fgFDMl1Wc08wSQlW778JLMCnATLwFHC3U7ujQ7zsA2Iq3lMotZlYe+v1eB+Y5534fp2XCJI5oiE9k/90HXI+35Hs3YD5eYhqO9wb8pHOuYl2jp/AS2al4dfyCwIfOueo+HJ6BV+Ovsrvw1nU6GW84chHe8ONQvMXxalrf52a8Ib/BeMntb865ThH8npOAf+EVHk4PbbszFONPgSMAw1vGo2Xo2tmbofb9gEvxitNeXnHCUOLdCJwYwfNLA6cEJbL/njCzaaFlFd7BW6vo1lAv4iG8nkUP51xPvB7RxWY2J/TmPA7IBE6r5tzD8IqkVnZP6BxfAtOBpWb2l9C6Q08AXUO9tap8ZmZPhOK7BW/kZOC+fsnQpJDdQLGZ5TnnUoAbgKvM7HMzW2Zm1+ENk44DWuIt6rjezFab2Sd4vbZPKp16aej3FKmRhvhE9l/48t6FwGozq6i6XLHuVhOge+hnc26vpcRS8XpVVc2y6whsiuA5w9uEP+fOKo7dM9nCzHaEYkmuot2+9Ao9x0ehYcQKTQFnZlucc/cCTznn7gDeA142szmVzrMZb10ikRopQYnsv8oL/lW3UnGjUNvBeEN74bZUc0yAqicmRPqcVSmpYlvFc1SOq6b3hIp9JwP5lfbtADCzW0PT98/CW4biY+fcn8zs/8LaJvF9UhWplob4RKInC6+n0szMss0sG29CwwNAdVOt84D2MYoPvOTVMuxxz0r7wxNYNt6KvB3Cfp+VwB3ACOdcN+fc34A1ZvaAmR0P/BlvEkm4dsRwZqTUXepBiUSJmZlz7m3gBefctXiTA+7Bm1yxrJrD5uHNuIuVOcCVzrn3gMbA3eydlAqAw5xz3c1sjXPur8CjzrliYAXeNamzgdvxfr+fADjnHgTS8HpblYf4BuKtSi1SI/WgRKLrMryp6G/ivVG3xKvMsK2a9u/izfaLldvwrlF9BbyA1xsKHz78B9AJyArN/LsJ+DfwHN5MwiHAGDNbaWa7gR/hTTFfiDc5YjHejccAOOcOxVvtubZuJJZ6TCvqisQR51wqsBo4zczm+xxOrQtNouhgZlf4HYvEP/WgROKImRXiXaO61u9YaltomvrFeOWTRPZJCUok/jwMDHSV5qbXA78Gnq+oOiGyLxriExGRuKQelIiIxCUlKBERiUtKUCIiEpeUoEREJC4pQYmISFz6f+UMiDtODncYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(results.T, label='milk')\n",
    "decorate(xlabel='Time (minutes)',\n",
    "         ylabel='Temperature (C)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
